Excitable media in open and closed chaotic flows 
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We investigate the response of an excitable medium to a localized perturbation in the presence of 
a two-dimensional smooth chaotic flow. Two distinct types of flows are numerically considered: open 
and closed. For both of them three distinct regimes are found, depending on the relative strengths 
of the stirring and the rate of the excitable reaction. In order to clarify and understand the role 
of the many competing mechanisms present, simplified models of the process are introduced. They 
are one-dimensional baker-map models for the flow and a one-dimensional approximation for the 
transverse profile of the filaments. 
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I. INTRODUCTION 



Excitable media [1-3] are extended systems exhibiting 
a variety of pattern formation phenomena. They are of- 
ten of chemical or biological nature, although they can 
also be found in other contexts [4,5] . Each spatial point 
in an excitable medium is described by a dynamical sys- 
tem in which activator and inhibitor variables can be 
identified. The activator displays some kind of autocat- 
alytic growth behavior, but the presence of the inhibitor 
controls it so that the dynamical system has a stable 
fixed point as unique global attractor. The essence of the 
excitability phenomenon is the presence of a threshold, 
such that if the system is perturbed above it, the system 
variables reach the stable fixed point only after a large 
excursion in phase space. This behavior usually appears 
when the activator has a temporal response much faster 
than the inhibitor, which then takes some time before 
stopping the growth of the activator. 

When different parts of a system are coupled diffu- 
sively, a local perturbation excites neighboring points, 
and as a result, the excitation propagates through the 
system as a wave (called autowave or front). This is a 
global phenomenon in the sense that all the points in 
the system will be reached by the wave and thus ex- 
perience the excitation-deexcitation cycle, but is non- 
coherent, since only a small part of the system (the front) 
is excited at each time. 

In many situations the excitable dynamics takes place 
in a fluid environment. One such example is the 
Belousov-Zhabotinsky (BZ) reaction [6,7], intensively in- 
vestigated in laboratory experiments. Another example 
is the population competition occurring in oceans or lakes 
between different plankton species: Truscott and Brind- 
ley [8] identified phytoplankton as the fast activator and 
zooplankton as the slow inhibitor in models of biological 



aquatic population dynamics. In such situations, differ- 
ent parts of the system interact not only via diffusion, 
but advective transport is present, and it can play also 
an important role [9]. One of the most obvious effects 
of stirring by the flow is that the concentrations would 
become more mixed and, for fast enough stirring, the 
whole extended system will behave homogeneously. An- 
other effect, in this case present for slow stirring, is that 
the fronts would be deformed by the flow and eventually 
be broken [10]. A study covering the full range of stir- 
ring intensities was performed in [11], in the framework 
of flows leading to chaotic advection [12]. Among other 
results, it was shown that there is an intermediate range 
of stirring for which the whole system excites coherently. 

In this Paper we further analyze the situation ad- 
dressed in [11], that is, excitable media under the effect 
of chaotic advection, and extend it by considering also 
stirring by open flows. As a striking result, we find sit- 
uations in which excitation persists indefinitely in the 
system when stirred by an open flow, whereas the excita- 
tion process is a transient both under closed flows and in 
the absence of stirring. Additionally, a number of simpli- 
fied one-dimensional models are introduced and used to 
gain insight and analytical predictions on the dynamical 
processes involved. We mention that studies in the same 
spirit than ours but for the different case of chemical re- 
actions of autocatalytic or bistable type can be found in 
[13]. 

The Paper is organized as follows: In Section II we 
present the basic framework and the chemical and two- 
dimensional flow models (closed and open) to be used. 
Section III describes numerical results for them. Section 
IV introduces one-dimensional simplified models that 
help to understand the above numerical results, and our 
Conclusions are presented in Section V. 
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II. REACTION-ADVECTION-DIFFUSION 
DYNAMICS 

A. General framework 

Let us consider N interacting species with concentra- 
tions Cj(r, t), i = 1, 2, ...N, transported by a flow v(r, t) 
that we assume incompressible. The governing reaction- 
advection-diffusion equations can be written as: 



da 

dt 



+ v • VCi = Ti (Ci, Cat; fc lf fc M ) + AV 2 ^, 

(1) 



where the Ti (Ci, Cjv; fci, /cm) describe the interac- 
tion dynamics of excitable type among the components, 
and the parameters fc; arc the reaction rates. Although in 
some realizations of excitable reactions the diffusion co- 
efficients can vary widely from one species to the other, 
as for example when it takes place in a gel medium, in 
liquid media diffusion coefficients are rather similar and, 
for simplicity, we take in the following Z), = D, Mi. It is 
convenient to adimcnsionalize Eq. (1) to have a clearer 
view of the processes involved. To this end, let us iden- 
tify typical scale, L, and speed U of the flow. A typical 
time scale is thus L/U, and we perform the change of 
variables: 



t 



t' = 



tu 



v 



V 

77' 



(2) 



We assume that the concentrations are already ex- 
pressed in some convenient dimcnsionlcss units, so that 
the reaction rates have units of inverse time, and use one 
of the reaction rates, say k\ = k, to define dimensionless 
reaction terms: 



T, 



Gi C\, Cat; ■ . ... 



ki km. 
fc'-' k 

k~ x Ti (Ci, Cat; k, k 2 , k M 



(3) 



With these changes, Eq. (1) reads 

dC 1 

+ v • VCi = Da G t (Ci, C N - e 2 , e M ) + 7^V 2 C 4 

(4) 

where the primes have been omitted for notational sim- 
plicity. We have scaled the reaction rates in terms of the 
first one ej = ki/k, i — 2, 3, M, and 



Da.f andPe.f 



(5) 



are the Damkohler and the Peclet number, respectively. 
The Damkohler number measures the reaction speed in 
terms of the advection, whereas Pe is the ratio of advec- 
tion to diffusion at scale L. The product PeDa measures 



the importance of reaction with respect to diffusion. We 
will be interested in the regime of large Pe, so that dif- 
fusion is negligible except at scales much smaller than 
system size, and explore a range of values of Da. We 
consider several two-dimensional models of flow, and also 
one-dimensional simplifications of them. Sensible com- 
parisons of the behavior under different flows would be 
facilitated by the introduction of the above adimensional 
framework, although perfect correspondence can not be 
expected when they are not dynamically similar. 



B. Reaction and flow models 

As a concrete example of reaction scheme of the ex- 
citable type, we focus on the FitzHugh-Nagumo (FN) 
model [1,2]. We note however that we expect all our qual- 
itative results to apply to the general class of excitable 
systems. In fact, some early results for closed flows [11] 
have been already checked for models of plankton dy- 
namics [14]. The FN model consists in a dynamics of the 
type (4) for two interacting species, of concentrations C\ 
and C 2l and reaction terms: 

Qi=f{Ci)-C 2 , f(C 1 )=C 1 (a-C 1 )(C 1 -l) (6) 
02 = e{d - 7 C 2 ) (7) 

£ (= £2) is the ratio between the two time scales. The 
FN model shows excitable behavior when e< 1 so that 
there is a separation between the fast evolution of the 
active component, C\, and the slow evolution of C2, the 
passive one or inhibitor. 

In a homogeneous system, Eq. (4) becomes C, = D&Qi . 
When C2 = 0, this dynamical system with Eq. (6) de- 
scribes dynamics of a bistable reaction, so that initial 
conditions for C\ below a threshold value a decay to- 
wards the unexcited or rest state C\ = 0, whereas initial 
conditions above a evolve to the excited state C\ = 1 in 
a time of the order of Da -1 . But Eq. (7) implies that, 
as soon as C\ grows above zero, the inhibitor C2 grows 
also (on a time scale a factor e slower) and as a result the 
excited state C\ s=s 1 is only a transient: in a time of the 
order of 



Cf(eDa) 



(8) 



Ci reaches C^f , which is the local maximum value of the 
function /(Ci), and then deexcitation occurs. After a 
time of the order of (eDa7) _1 , during which the system 
can not be excited (the refractory state) C\ and Ci return 
back to the fixed point or equilibrium value C\ = Ci = 0. 
In the following we will use the values e = 1CP 3 , 7 = 3.0, 
and a = 0.25 for the parameters in the local FitzHugh- 
Nagumo dynamics. In this case, C^ 1 w 0.1. 

The flow is assumed to be imposed externally so that 
the chemical dynamics has no influence on the velocity 
field, v(r, t). We consider two different kinds of chaotic 
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flows which, in other contexts, are known to behave 
rather differently: closed and open flows. In the first 
situation, fluid particles remain in a bounded region of 
space, and the flow produces mixing in the whole fluid. 
In the case of open flows, fluid particles enter the system 
and typically, after some time, they leave it. Interesting 
situations arise when, as an effect of the stirring, there 
are special orbits never leaving the system. For hyper- 
bolic chaotic flows, such orbits form a fractal set of zero 
measure, the chaotic saddle [15] with stable and unsta- 
ble manifolds. When a set of particles is released on 
that flow, most of them leave the system rapidly. But 
those on trajectories coming close to the stable manifold 
of the chaotic saddle become attracted by it and remain 
longer in the system. They finally escape the saddle, at 
a characteristic escape rate k, tracing closely its unstable 
manifold. 

As a simple two-dimensional and incompressible veloc- 
ity field, we consider an archetype of closed chaotic flows, 
the motion generated by two alternating sinusoidal shear 
flows oriented along the x and y direction for the first 
and second half of a flow period, T, respectively [16]: 

x A^/T , \ . (2iry , \ 

v x (x,y,t) = —0^--tmod Tj sin I — + fa I 

A / T\ ( 2nx \ 

v y (x,y,t) = —e[tmodT- -j sin ( — + fa+i j, (9) 

where is the Heaviside step function. Note that all 
the geometric details of the stirring by the flow depend 
on the parameter A, while T sets the speed of stirring 
without altering the trajectories of the fluid elements. 
In order to avoid KAM tori acting as transport barri- 
ers, typically present in time-periodic flows, the velocity 
field is made aperiodic by introducing a random phase 
cj>i, which takes on independent values in each half pe- 
riod, and is uniformly distributed in the range [0, 2ir]. 
The fluid is confined in a square of lateral size L with 
periodic boundary conditions, so that the flow is closed. 
L and T fix space and time scales for the flow, and L/T 
is a typical velocity. Thus we can adimensionalize Eq. 
(1) in terms of these quantities, so that Da = kT and 
Pc = L 2 {DT)~ l . This leads to Eq. (4), written in units 
for which L = T = 1. We set the remaining adimcn- 
sional parameter A/L = 0.7, for which the flow is nearly 
ergodic. The numerically computed Lyapunov exponent 
is fi w 1.66/T. 

As a simple example of open flow, we take a blinking 
vortex-sink system [17,18] consisting of two alternately 
opened point sinks in an unbounded two-dimensional do- 
main. Around each sink, the velocity field is a combi- 
nation of a point-vortex and a point-sink given by the 
complex potential 

w(z) = -(Q + iK)ln\z- z s \ . (10) 

z s gives the position of the sink in the complex plane 



{z, z G C}, so that z — z s — re 1 ^ defines polar coordi- 
nates (r, </>) around the sink. The imaginary part of w(z), 
^ = — Klnr — Q(j> is the streamfunction, from which the 
fluid particle equations of motion are 

• - - 9. 

r def) r 

d$> K 

r4> = -— = -. 11 

or r 

The fluid trajectories can be explicitly integrated: 
r(t) = y/r$-2Qt 

0(t)=to-|]n^. (12) 

The fluid particles come from infinity following logarith- 
mic spirals of circulation given by 2nK. The flow is 
incompressible everywhere, except at the sink core z s , 
where an area of fluid 2ttQ disappears per unit of time 
(the trajectories in a circular region of this area around 
z s have their trajectories undefined after one time unit 
because of (12) and they should be understood to leave 
the system): we have thus an open flow. The motion is 
not chaotic if z s is a unique static point. The blinking 
vortex-sink flow consists in considering the active sink 
position to be at z s = (0,d) during half a period T/2, 
and at z s = (0, —d) during another half a period. This 
corresponds physically to opening and closing alterna- 
tively two sinks separated by a distance 2d. Typical 
space, time, and velocity scales are d, T, and d/T, re- 
spectively. From them, Da = kT, and Pe = d 2 (DT)~ 1 , 
and we can use Eq. (4) in the units in which d = T = 1 . 
The flow is fully characterized by the adimensional sink 
strength 77 = QT/d 2 and the ratio of vortex to sink 
strength £ = K/Q. We take 77 = 1 and £ = 10. For 
these parameter values, the Lyapunov exponent on the 
saddle and the escape rate from it are [18] \i = 2.19/T 
and k — 0.54/T, respectively. 

III. NUMERICAL RESULTS 

The numerical integration of the reaction-advection- 
diffusion problem has been carried out on a square grid 
of 1000 x 1000 points, with grid size Ax, by using a sim- 
ple semi-Lagrangian scheme for the transport processes 
combined with a fourth-order Runge-Kutta method for 
the time integration of the local chemical dynamics. The 
semi-Lagrangian advection step at time t consists in 
calculating, from any gridpoint, a time-backwards La- 
grangian trajectory for a time At. Then the concen- 
trations at this fluid element are calculated by bilinear 
interpolation from the concentrations on the grid at time 
t — At, and these concentration values are then assigned 
to the starting gridpoint at time t. The interpolation step 



3 



introduces an effective diffusion D w (Ax) 2 /At, which 
limits the maximum Pe number we can attain. Since the 
numerical diffusion is not uniform in space, we also in- 
clude an explicit diffusion step corresponding to the same 
Pe number. 

Initially the system is in the homogeneous steady state, 
C\ = C2 = 0.0, and then it is perturbed by a localized 
Gaussian pulse in the concentration of the activator com- 
ponent 



d(x, y,t = 0) = C cxp (-(.x 2 + y 2 )/2l 2 ), 



(13) 




where Co is chosen to be larger than the excitation 
threshold a — 0.25, and the size of the perturbation, lo, 
is much smaller than the system size, L = 1. The depen- 
dence of the results on the particular values of Co and 
Iq will be discussed later. The inhibitor component, C%, 
is not perturbed initially. We study the response of the 
system for different values of the adimensional reaction 
rate Da, keeping the rest of the parameters fixed. In the 
absence of flow (Da = 00) the initial condition (13) pro- 
duces a circular ring of excitation (a target wave; in fact 
the structure of the target wave is such that the excita- 
tion ring is followed by a refractory ring) that expands 
in radius until reaching the system boundaries. We will 
see that this behavior is strongly modified at finite Da. 




FIG. 1. Activator concentration at different times under 
stirring by the closed flow, for Da = 300 and Pe » 1000. The 
initial condition had Go = 0.5 and lo = 0.01. The sequence 
runs from top to bottom, and then from left to right. The 
total time lapse is 3 periods of the flow. Dark and clear gray 
level indicates, respectively, low and high activator values. A 
noncoherent process of global excitation is seen. 



A. Closed flow 



The model (9) is integrated numerically on the unit 
square, L = 1, with periodic boundary conditions, from 
the initial condition described above. The resolution is 
thus Ax = 10~ 3 . We use At = 10~ 3 (in units of the flow 
period) and thus Pe « 1000. 



Snapshots of the spatial structure for Da = 300 and 
Da = 25 are shown in Figs. 1 and 2, respectively. In 
the first case, the localized perturbation gives rise to an 
excited patch (with its interior in the refractory state) 
that is the stirred version of the circular wave front that 
would be produced in the absence of flow. The patch 
is elongated into a convoluted filamental structure by 
the chaotic flow and eventually visit all the points of the 
system. The filaments have a characteristic double-line 
structure with refractory area in the center. The exci- 
tation is global, in the sense that all the points of the 
system have become excited at some moment, but is not 
coherent, since only a part of the system is excited at a 
given moment, being the rest in the refractory or in the 
equilibrium state. 
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pie, if Co = 1, and Iq = 0.1, 0.05, and 0.02, one finds 
Da c w 12.0, 14.2, and 16.3, respectively. 




FIG. 2. Same as Fig. (1), but at Da = 25. The interval 
between snapshots is of 1.2 time units. A process of coherent 
global excitation is seen. At longer times (not shown) the 
system returns homogeneously to the unexcited state. 

In a range of smaller Da numbers, a qualitatively dif- 
ferent phenomenon occurs: The initial patch is again 
stretched into a growing filament, but now the filaments 
are thinner that prevents the formation of refractory re- 
gion within them. This results in a coherent global exci- 
tation when the filaments fill up the whole system. Fig. 
2 is a representative example of this situation. Once fully 
excited, the system remains homogeneous and its subse- 
quent decay to the unexcited state occurs everywhere at 
the same time. In this second part of the dynamics, mix- 
ing becomes irrelevant since there is nothing to mix in an 
homogeneous configuration. 

At still smaller Da (faster stirring or slower chemistry) , 
a sharp transition to a new dynamic regime occurs: be- 
low a critical Da number, Da c the excitation dies with- 
out propagating significantly: dilution is fast and dom- 
inates over the growth rate of the activator. By dilu- 
tion we mean the process of mixing of the excited patch 
with the surrounding unexcited fluid, leading to the de- 
creasing and elimination of excitation in the patch. This 
mixing is originated by the diffusive flux from inside the 
patch (high value of the activator) to outside (low ac- 
tivator value). The value of Da c depends only slightly 
on the details of the initial condition (as long as Co re- 
mains suprathreshold and Iq much smaller than system 
size and above some diffusion-controlled minimum size, 
of the order of Eq. (22), discussed below). For exam- 
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FIG. 3. Time evolution of the mean concentration at var- 
ious values of Da. Co = 0.5 and lo — 0.01. Time is in units 
of T. By further decreasing Da, this time evolution changes 
suddenly, at Da c , to a fast decay of the initial condition that is 
indistinguishable, at the scale of this plot, from the horizontal 
axis at zero mean concentration. 
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FIG. 4. The maximum value attained by < C\ > as a func- 
tion of Da. In all cases the initial condition had Co = 1, and 
lo was as indicated. 

Figs. 3 and 4 summarize the different situations. Fig- 
ure 3 shows the time dependence of (Ci), the space- 
average of the concentration of Ci, and Fig. 4, the maxi- 
mum value attained by (Ci). For increasing Da the co- 
herence of the global excitation is gradually lost so that 
the maximum value of the average concentration is below 
the one corresponding to the fully excited state. 



B. Open flow 

Since an infinite domain can not be simulated easily in 
the computer, a square domain of size L/d = 6 is consid- 
ered instead. With a lateral discretization of 1000 points, 
this leads to Ax = 6 x 10" 3 . We use At = 2 x 10~ 3 (in 
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units of the flow period) and thus Pe « 56. Concentra- 
tions at the boundaries are kept at the fixed-point values 
C\ = C-2 = 0. The interior of the domain is initial- 
ized also in this state except for the perturbation in C\, 
Eq. (13), located in the middle position between the two 
sinks. 

For small Da, as in the closed flow case, the per- 
turbation is diluted by the flow before significant wave 
propagation occurs, and the excited material leaves soon 
the system through one of the sinks. By increasing 
Da, a sharp transition to a new regime occurs: In re- 
sponse to the persistent arrival of unexcited reactants 
from the boundaries, and despite the continuous loss of 
fluid through the sink, a steady pattern of excitation is 
permanently sustained by the autocatalytic behavior of 
the activator. An example of the time evolution is shown 
in Fig. 5, and a snapshot in Fig. 6. The excited pattern 
closely traces a fattened version of the unstable manifold 
associated to the chaotic saddle of this open flow [18], 
and as this manifold, it fluctuates periodically in time. 
The value Da c « 14.5 at which this transition occurs is 
in the range of the values obtained for the closed flow, 
thus suggesting that the mechanism for it is rather local 
and independent of the details of the flow. 



If 




^^^^^^^^^^ 

















FIG. 5. Time evolution of the activator concentration un- 
der the open flow. The snapshots are shown every 0.7 time 
units, from top to bottom and then from left to right, start- 
ing 0.7 time units after the initial perturbation, Da = 50 and 
Pe m 56. After the last time shown, excitation is maintained 
indefinitely in the system, and the long-time pattern, which 
follows the shape of the unstable manifold of the chaotic sad- 
dle, repeats periodically in synchrony with the flow. 




FIG. 6. snapshot of the activator concentration maintained 
in the system at long times, for Da = 50. 



The filaments in the excited pattern fatten up with in- 
creasing Da. Suddenly, a new regime is reached above a 
second critical Da ss 90.5: the excitation initially accu- 
mulates at the unstable manifold of the chaotic saddle, as 
before, but this is just a transient that is followed by an 
irregular recovery of the equilibrium (C\ = C2 = 0) state 
everywhere. This new regime has some analogies with 
the large Da behavior under the closed flow, for which a 
noncoherent excitation occurred (and the transition value 
of Da is of the same order), but here it appears suddenly 
as a function of Da, and the excitation does not visit the 
full system but remains close to the unstable manifold of 
the chaotic saddle. 

The whole behavior is summarized in Figs. 7 and 8, 
where the time evolution of the mean value of the acti- 
vator, and its asymptotic long time value (not the maxi- 
mum value as in Fig. 4), is plotted versus Da. There is 
a range of Da in which a finite amount of excited fluid 
remains permanently in the system, despite the openness 
of the flow. 

The existence of critical values of Da or equivalently, 
of critical stirring rates can be seen to be a consequence 
of the competition between a number of processes. In the 
closed flow case, advection and diffusion tends to homog- 
enize and dilute the excited patch, while the excitable 
dynamics increases the local concentration of the active 
component wherever the excitation threshold is exceeded. 
In the open flow, there is an additional factor: the escape 
rate of fluid particles from the system. In order to gain 
further insight we attempt to separate the essential in- 
gredients that contribute to the observed behaviour, and 
consider reduced models of the problem. 
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FIG. 7. Time evolution of the average concentration 
< Ci > for three values of Da. Co = 0.5, lo = 0.05. 



Even for arbitrary initial conditions the concentrations 
are rapidly homogenized along the y direction by the 
repeated stretching, and after a short time the one- 
dimensional description becomes relevant. This may be 
regarded as a general feature of transport problems in the 
presence of stirring, since one can associate local stretch- 
ing directions to any point of the flow and the problem 
can be reduced to the description of the fhamental struc- 
ture in the transverse direction. In the one-dimensional 
formulation the baker transformation T acts by replac- 
ing the concentration field by two copies compressed by 
a factor of two placed next to each other. To better rep- 
resent the process of filament folding, the left half is not 
a copy but the mirror image of the right half. 
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FIG. 8. Long-time asymptotic value of the average concen- 
tration < Ci > as a function of Da. It is non-zero in a range 
of intermediate values. 



T: x^>Tx = x/2,{2-x)/2. 

d{x) - Td(x) = C^r^x) (14) 



Numerical simulations on the unit interval with periodic 
boundary conditions of the one-dimensional FitzHugh- 
Nagumo system with the baker transformation applied 
at discrete times (t = nT, n = 1, 2, ..), and diffusion and 
chemistry acting between them, show qualitatively simi- 
lar regimes to the two-dimensional closed flow presented 
above, including the transition to global excitation that 
occurs in an intermediate range of Da = kT. 



IV. REDUCED MODELS 



A. One-dimensional baker model for the closed flow 



The main effect of chaotic advection is to stretch and 
fold fluid elements producing the filamentary patterns 
visible in the figures of the last Section. Perhaps the 
simplest model of this is the so called baker transforma- 
tion. A single action of the baker transformation on the 
unit square can be described as a stretching along the y 
axis by a factor of two, followed by compression by a fac- 
tor of two along the x axis. Then the resulting rectangle 
is cut into two pieces of unit length along the y direc- 
tion and placed back on the unit square. This model 
of chaotic advection neglects spatial non-uniformities of 
the stretching and curvature of the filaments, present in 
a general flow. Nevertheless, since it is a discrete-time 
map, it is strongly non-uniform in time. 

If the initial perturbation is taken to be homogeneous 
in the y direction this is preserved by the baker trans- 
formation and the problem becomes one-dimensional. 
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FIG. 9. Spatiotemporal evolution of Ci for Da = 1 and 
Pe = 1000 (upper left), Da = 10 and Pe = 1000 (upper 
right), Da = 40 and Pe = 1000 (lower left), and Da = 50, 
Pe = 400 (lower right) under the baker model. Space is in 
the horizontal direction, and time runs in the vertical from 
bottom to top. Darker grey represents smaller values of C\. 
The discontinuities appear at each application of the baker 
map, i.e., at times T, 2T, 3T, etc. 

Time evolutions of the activator spatial structure are 
shown in Fig. 9, for three different values of Da, and two 
of Pe = L 2 /DT. The excited regions can be interpreted 
as transverse cuts through filaments. The number of fila- 
ments is doubled by each action of the baker map, while 
the decreasing of their width may be, or may be not, com- 
pensated by the effect of excitable growth. When growth 
is slow (upper left panel in Fig. 9), the filaments become 
narrower until a point in which diffusive mixing with the 
surrounding unexcited fluid destroys them and excitation 
disappears. By increasing Da (upper right), the filament 
width reaches a minimum nonvanishing value with its 
central C\ concentration value well above the threshold 
a. The effect of the baker map is here to join together a 
number of filaments until diffusion homogenises the dis- 
tribution. Since the homogenized value of C\ is above 
the threshold, a coherent excitation follows. After some 
time the excitation disappears homogeneously and the 
system returns to the nonexcited state. 

In the case of large Da (lower panels), the reaction 
is fast enough to approach the refractory state in the 
middle of the filaments before significant compression. 
The filaments thus acquire the double-hump structure 



that was also seen in the full two-dimensional simula- 
tion (Fig. 1). This fact works against the possibility of 
a coherent excitation, and there are two mechanisms by 
which the noncoherent excitation dies. At large enough 
Pe (lower left panel of Fig. 9), the excited parts of the 
filaments are narrow, and the periodic contraction pro- 
duced by the baker map eventually bring them below a 
width such that diffusion can eliminate them, in a way 
similar to what happens with the full filament at small 
Da (but here there is around abundant refractory mate- 
rial that helps the process). When Pe is decreased, the 
excited parts of the filaments travel faster, and collisions 
leading to filament annihilation (because of the refrac- 
tory material arriving after the excited filament) are the 
main mechanism killing the excitation (lower right panel 
of Fig. 9). 

The phenomenology found here is fully consistent with 
the numerical simulations of Section III. But now, in 
addition to having a much simpler numerics, the mecha- 
nisms are easier to identify. Thus, stretching and folding, 
the characteristics implemented in the baker map, are 
enough to understand the effects of chaotic advection on 
excitable advection-reaction system. We stress however 
that the Da values at which the different transitions oc- 
cur are of the same order, but not identical, to the ones 
found in the two-dimensional models. This was expected 
since there is no complete dynamical similarity between 
the present baker model and the flow models of Sect. III. 

In the baker map model the spatial structure is, by 
construction, periodic with period L/n,n = [t/T]. Thus 
the evolution of the system can be fully described by 
solving the same problem on an interval compressed by 
a factor of two at times t = nT with periodic boundary 
conditions. This suggests, in general, that the evolu- 
tion of the system can be captured by focusing on the 
transverse profile of a single filament subject to a typical 
stretching, and taking into account the decreasing sep- 
aration between the filaments by appropriate boundary 
conditions. In fact, the main mechanism controlling the 
final homogenized value is the competition between the 
compression by the flow, and the tendency to expansion 
due to reaction-diffusion. This mechanism is better an- 
alyzed by considering an isolated single filament, as will 
be done in Subsection IV C. 



B. One-dimensional baker model for the open flow 

As in the closed flow case, we can implement the es- 
sentials of chaotic advection in open flows: stretching, 
folding, and escape, by a one-dimensional version of the 
open baker map. At times nT, n = 1,2, the unit in- 
terval (L = 1) is compressed a factor of three; two copies 
of the resulting compressed configuration are placed back 
into the initial square (one of them with orientation re- 
versed) and the remaining third is filled with unexcited 
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material (Ci = C2 = 0). This represents the loss of one 
third of the fluid per map step, and its substitution by 
fresh reactants. Standard diffusion with periodic bound- 
ary conditions, and FitzHugh-Nagumo dynamics act be- 
tween successive applications of the map. 




FIG. 10. Concentration of C\ under the open baker map 
for Da = 30 (left) and Da = 50 (right). Space is in the 
horizontal direction, and time in the vertical, running from 
bottom to top. Darker grey levels correspond to lower values 
of Ci. Pe= 10 4 . 

The phenomenology observed is again qualitatively 
consistent with the two-dimensional simulations. For 
small Da the initial excitation is diluted before signifi- 
cant propagation. At larger Da (Fig. 10, left panel ), 
the excitation approaches the chaotic saddle of this map, 
which is a standard Cantor set, and covers it with a fi- 
nite width. A dynamic equilibrium is reached between 
filament merging and filament replication, so that the 
excitation is maintained indefinitely in the system. In- 
creasing further Da leads to a second transition to a situ- 
ation in which the excitation finally disappears, in much 
the same way as in the closed flow. Again this happens 
when the filaments begin to develop the refractory state 
in its interior, and it may occur by two mechanisms: the 
one shown in the right panel of Fig. 10 which involves 
complex filament interaction, or simply the repeated con- 
traction of the narrow excited parts at both sides of the 
refractory center. This last mechanism dominates at very 
large Pe. 

Both in the open and in the closed flow case, the mech- 
anisms leading to transitions and qualitative changes in 
the excitation behavior seem to be linked to properties 
of individual filaments, namely, the existence of a mini- 
mum width below which the filament disappears, and the 
development of a double-hump shape by increasing Da. 
Both phenomena can be understood to great detail by 
focusing on the behavior of an isolated filament. 



C. A one-dimensional filament model 

We can address the analysis of the one-filament prob- 
lem by replacing the flow by a time-continuous stretching 
in a pure strain flow, v x = —Ax, v y = Ay. This has in 
common with chaotic advection and with the baker map 
the local contraction and expansion along special direc- 
tions, although it misses completely the folding behavior 
that leads to filament interaction at long times. Accord- 
ing to the discussion above the relevant dynamics is along 
the convergent direction (x) of the flow. Thus, we pro- 
pose that the evolution of concentrations Ci(x,i) in a 
chaotically advected excitable medium can be described 
locally by 

®c i -\x-?-C i =F i (C u ...,C N )+D^Lc h (15) 
at ox dx 

where A is the strain due to the flow [19,11,13]. In gen- 
eral, the strength and direction of the stretching fluctu- 
ates in space and time. Thus a suitable prescription for 
fixing a unique A should be established. This issue will be 
discussed later. A way to take into account multifilament 
situations in the framework of Eq. (15), is to impose pe- 
riodic boundary conditions on an exponentially shrinking 
interval of length L = cxp (—At), taking into account the 
decreasing interfilamental distance. But we will consider 
here the case of an isolated filament in a large (ideally 
infinite) one-dimensional domain. 

We note that the one-dimensional model (15) does not 
conserve the amount of fluid on the line. This is more 
clearly seen be rewriting it as 

l Ci + T x {- XxCi ~ D lk c )= (16) 

FiiCi, Cjv) — ACj . 

Whereas the left hand side is clearly written in a flux- 
conservative form, the term — XCi in the right hand side 
represents fluid escape from the line at a rate A. The 
reason for that is that Eq. (15) comes from a strain flow 
in which there is motion along the y axis, and the term 
XCi is simply the flux in that direction for concentrations 
homogeneous along y: AC,(x) = d y (\yCi(x)). 

The pure strain flow has a time scale, A" 1 , that can 
be used to adimensionalize times, but there is no typical 
length scale. However, we can measure lengths in units 
of the diffusion length y/D/X and then Eq. (15) becomes 

- x JrC, = Da&(Ci, ...,C N ) + (17) 
at dx dx 

with Da = k/X, t = Xt, and x = x{D/X)~ 1/2 . Thus, we 
can always set Pe = L 2 X/D = 1 by choosing the units 
of L or, in other words, the only effect of variations of 
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diffusion strength in this model is a change in spatial 
scale. Qualitative changes can only occur by varying Da, 
that is, by changing k or A. This is clearly a limitation 
of the model and tell us that it can only be trusted in 
regions where there are well separated filaments of size 
much smaller than characteristic spatial scales of the ve- 
locity field (which arc neglected when assuming a pure 
strain). We expect this to be a reasonable global ap- 
proximation at sufficiently large Pe. Other phenomena 
neglected by this one-dimensional model are strain inho- 
mogeneities and departure from one-dimensionality. 

In this Section we mainly present our results in terms 
of the parameter Da of Eq. (17), but eventually we would 
need to return to the units of Eq. (15), where the indi- 
vidual processes and scales are more easily identified. In 
the search of clarity, quantities representing lengths will 
be marked with an overbar when measured in the units 
of Eq. (17), that is, in units of the diffusion length. We 
note also that changing the strain A in Eq. (15) changes 
Da = k/X and also the units of space and time in Eq. 
(17). 

Numerical solution of Eq. (17) for Da not too small 
reveals that its long-time attractors are steady pulses of 
excitation concentrated near the origin. They can be 
interpreted as transverse cuts of the filaments observed 
in the two-dimensional models. Examples are shown in 
Figs. 11a) and lib), where we plot the C\ and C 2 con- 
centration fields, respectively, for different values of Da 
(the insets will be discussed below). The steady finite 
width of the filaments arises from compensation between 
the contracting tendency of the strain, and the expand- 
ing tendency of the combined effect of diffusion and re- 
action. A simple quantitative argument [11] formalizing 
this consists in identifying the equilibrium half-width of 
the filament solutions of (15), w s , as the distance to the 
center at which the strain speed Xw s exactly compensates 
the speed the front would have in the absence of strain, 
Vf. Since this velocity may be calculated for small e [2]: 
Vf = (1 — 2a)yJ Dk/2, w s is obtained from Xw s = Vf. In 
the adimensional space units and parameters of Eq. (17) 
it reads: 

w s = (l-2a)y^. (18) 

The existence, for Da above (or A below) a given value, 
of these steady filaments with finite width provides an ex- 
planation for most of the phenomenology discussed in the 
previous sections. In two-dimensional situations, the fila- 
ment will maintain its transverse shape while it expands 
in the longitudinal direction. After repeated folding, it 
will cover the whole system in the closed flow case, or 
cover the unstable manifold of the chaotic saddle in the 
open flow case. What will happen latter will depend on 
the interactions between different parts of the excited fil- 
ament, or on its response to strain fluctuations. 



a) b) 




FIG. 11. a) Ci profiles for stable one-humped solutions of 
Eq. (17). In the inset the corresponding unstable solutions. 
Solid line is for Da = 12.99, close to the disappearance of the 
filament solution, dotted-line for Da = 16.77 and dashed-line 
for Da = 25. b) The same as a) but for the C2 concentration 
field. All the curves are symmetric with respect to x = 0, so 
that only positive x values are plotted. 



We observe that thejsteady-filament stable solution 
disappears for Da < Da c w 12.5). This provides an 
explanation for the absence of excitation in the two- 
dimensional simulations below a critical Da: as far as 
the results of the one-dimensional model can be extrapo- 
lated there, a growing filament state can not be reached 
at small Da because a steady (non-decaying) solution of 
the filament profile docs not exist. To better understand 
the disappearance of the filament solution, we note that, 
in addition to direct numerical simulation, an alterna- 
tive way of finding the steady filaments is to solve by 
a shooting method [20] the steady state version of Eq. 
(17) which is obtained by setting c%Cj = 0. With this 
method one can obtain all the steady solutions, not only 
those that are dynamically stable. It turns out that, in 
addition to the excited filament (and to the trivial homo- 
geneous solution Ci(x) = C2(x) = 0) , there is another 
pulse solution, which is dynamically unstable. J?his un- 
stable solution is shown, for several values of Da, in the 
insets of Fig. 1 1 . It contains a marginal amount of excita- 
tion, in the sense that initial conditions with slightly less 
excited material evolve towards the stable homogeneous 
state, and initial conditions slightly more excited lead to 
the stable excited filament. It represents the unstable 
point in function space at which the activator growth ex- 
actly compensates the diffusive flux towards the exterior. 
In Fig. 12 we plot the width of the stable and unstable 
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steady filament solutions of Eq. (17). Here it is clear 
that the disappearance of the stable filament arises from 
collision with the unstable pulse in a saddle-node bifur- 
cation. Physically, increasing strain reduces the width 
of the stable filament, so that it approaches the unsta- 
ble one, which is the limit below which excitation decays. 
The saddle-node bifurcation will occur when both widths 
are equal. 



CUx) = c + 



C+-C- 



I _ tanh 2 (^-) 



c+ 



with 



C± = ^(l + a)±W^(l + a) 2 -2a. 



(20) 



(21) 



C- is the maximum concentration at the center of the 
pulse, and its half- width is given by w u = 2y/D/(ak), or 
in the adimensional units of Eq. (17): 
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FIG. 12. Half-widths measured at the level of the excita- 
tion threshold (Ci = 0.25) for the Ci stable and unstable 
filament solutions of Eq. (17) in terms of Da (log-log scale). 
Circles show the numerical values and the lines the analytical 
curves: solid-line from Eq. (22) and dashed-line from Eq. (18). 
The inset shows the same numerical widths but in the dimen- 
sional units of Eq. (15) (we use k = 1 and D — 10~ 5 ) as 
a function of A, to stress the insensivity of the width of the 
lower (unstable) branch to the strain A. 

Since the width of the unstable pulse around x = 
is rather small, we expect strain effects to be of minor 
importance in determining its shape, at least when Da is 
not too close to Da c . This is confirmed by the inset in 
Fig. 12, and also in the inset of Fig. 13 (to be discussed 
later). Hence, we can analytically estimate the shape of 
the unstable pulse and Da c in the following way: Since 
the amount of inhibitor is small everywhere for this so- 
lution (see the inset in Fig. lib)), and since we are in- 
terested in the situation e — > 0, we can approximate Eq. 
(15) (for A w 0) by 



fcCi(a-Ci)(Ci 



The unstable pulse, Cf(x), is the solution homoclinic to 
C\ = 0. After multiplication of (19) by d x C\{x), inte- 
grations with respect to x, and application of the proper 
boundary conditions, one finds [21] 



v 7 ; 



aDa 



(22) 



In the inset of Fig. 13 we compare the analytical curve 
from Eq. (20) with the numerical values. We see that 
they are similar but not identical. We have checked that 
the reason for the discrepancy is the finite value of e 
(e = 10~ 3 ). We have also performed calculations with 
smaller values of e and seen that for the values of Da 
in the figure, the approximate analytical solution (20) 
and the numerical one become virtually identical when 
e < 10~ 5 . 




0.01 0.02 0.03 0.04 

x(D/X)" 2 



FIG. 13. Two-humped filament solutions for the C\ field. 
Solid-line is for Da = 45.45, when the two humps begin 
to develop. Dotted-line Da = 50, dashed-line Da = 83.33, 
dashed-dotted line Da = 100 and dashed-double-dotted line 
for Da = 125. In the inset we show the unstable solutions, for 
Da = 45.45, 50, 52.63, and 55.55. In the scale of the plot, all 
of them collapse into the same solid-line curve, when plotted 
in terms of the physical space units of Eq. (15), as done here 
(we use k = 1 and D — 10~ 5 ). The dashed-line corresponds 
to the analytical solution given by Eq. (20). 



11 



In any case, since the widths of analytical and numeri- 
cal unstable pulses are very similar, we can estimate Da c 
by equating the above expressions, (18) and (22), found 
for them: w s = w u , with the result 



Da c = 



2V2 



(1 - 2a)y/a 



(23) 



For a = 0.25 this gives Da c w 11.31, tliat compares well 
with the numerically obtained value Da c w 12.5 (See Fig. 
(12)). 

The saddle-node disappearance of the filament solu- 
tions in this one-dimensional filament model clearly gives 
an explanation for the sudden disappearance of exci- 
tation propagation at small Da in the two-dimensional 
models discussed in Sect. Ill, and in baker-like models. 
To make the connection more quantitative, within the 
uncertainty given by the observed weak dependence of 
Da c on characteristics of the initial perturbation (Fig. 
(4)), one needs to identify the effective strain A in Eq. 
(15). If one identifies A w T _1 , which is a reasonable 
measure of the__strain in the models of Sect. Ill, then 
one has Da = Da and finds good quantitative agreement 
between the critical values of the Damkohlcr number for 
the filament model and the full two-dimensional simula- 
tions both in the open and in the closed flow case. But 
it should be said that, since the filaments are being ad- 
verted by the flow, a more consistent choice for A would 
be the Lagrangian mean strain, given by the Lyapunov 
exponent /j, of the advection dynamics. In the open flow 
case, the Lyapunov exponent on the chaotic saddle would 
be the analogous choice. These elections have been shown 
to be quantitatively successful in other situations [13]. 
With the values of /j, stated in Sect. II B, this leads to 
Da = 1.66Da for the closed flow and Da = 2.19Da in 
the open case. Now the agreement has deteriorated. The 
effect of strain inhomogeneities may be rather important 
when the filaments are wide and have some diffusive mo- 
tion, since then they can feel effective strains different 
from the Lagrangian one corresponding to a fluid particle 
at its center. Other effects related to the reduced dimen- 
sionality are also at play, since quantitative departures 
from the two-dimensional simulations appear already for 
the baker-like models of Sects. IV A and IV B. Thus, 
one concludes that the one-dimensional filament model 
needs to be improved to provide systematic quantitative 
predictions on the behavior of reacting systems, but it 
does a very good job in identifying the basic mechanisms 
and qualitatively modeling them. 

Still remaining to be discussed within the framework 
of this Section are the qualitative changes of behavior oc- 
curring in the two-dimensional simulations at large Da: 
the progressive loss of coherence in the closed flow case, 
and the sudden disappearance of the persistent pattern 
in the open flow case. These phenomena were more or 
less coincident with the appearance of a double hump 



structure in the filaments. Fig. 13 shows that indeed the 
stable filament solutions of the one-dimensional model 
(17) develop a double- humped structure for Da > 45. 
One can understand this by noticing that the front so- 
lutions of Eq. (15) for A w have a finite width lim- 
ited by the time during which excitation persists in the 
fluid particles (8), i.e., r e = C^ f (eDa) _1 , in units of A -1 . 
The width of the front is given in first approximation 
by Wf ~ VfT e . Interaction with the back of the front 
changes C^ 1 to a smaller value C/, which is the solution 
of an algebraic equation [2]. For a = 0.25, C/ w 0.067. It 
is reasonable to expect that, when the total width of the 
filament, 2w s , exceeds twice the width of the front 2wf, 
the filament will become unexcited in the middle. This 
argument would need corrections by the strain influence 
on Wf and by the fact that the strain velocity in the mid- 
dle of the filament is smaller than in the front (this would 
imply a shorter time of excitation, or smaller C/). But 
in any case, this simple argument gives for the transition 
to two- humped filaments the condition w s w Wf, that 
is Da w C/e -1 w 67, to be compared with the actual 
numerical value Da w 45. 



We mention that, the transition from the unimodal 
steady solution to the two-humped steady solution in the 
range 50 < Da < 70 is associated with a complex bifurca- 
tion scenario in which different stable solutions coexist. 
For Da > 70, additional asymmetric steady stable solu- 
tions to Eq. (17) appear. These are similar to the pulse 
solutions obtained without strain, but stopped by the 
flow. They have an excited head, and a refractory tail. 
An example is shown in Fig. 14. The symmetric two- 
humped filaments found before, which remain linearly 
stable, can be thought as bound states of two asymmet- 
ric ones. 
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FIG. 14. Solid line is an asymmetric steady solution of Eq. 
(17) for Da = 200. The higher curve is C\ and the lower 
one C2, multiplied by a factor of 5 to be visible in this scale. 
The dashed lines are symmetric solutions similar to the ones 
in Fig. 13, for Da = 200 (upper curve, Ci, lower curve, C2 
multiplied by 5). Numerically we find that both solutions are 
simultaneously stable, each with its own basin of attraction. 

The simple model studied in this Section has allowed us 
to qualitatively understand the individual filaments seen 
in the full two-dimensional simulations, and in the baker 
model, to a great detail. What is completely missed here 
is the interaction between different filaments, or parts of 
the same filament. We have studied filament collisions 
in the context of Eq. (17) by initializing it with differ- 
ent combinations of displaced filament solutions. The 
analogy with the collisions in the models of the pre- 
vious sections is far from complete, since here all the 
filaments evolve in the same simple velocity field —Ax, 
whereas in real multifilament situations, each filament 
has been created around its own local strain. Neverthe- 
less we have observed that collision between symmetric 
one-humped filaments leaves at long times a single cen- 
tered one-humped filament, and collisions between two- 
humped filaments annihilates half of the humps, leading 
again to a single two-humped filament as the final state. 
The asymmetric front-like filaments annihilate when col- 
liding front to front, and bind in a two- humped filament 
when colliding tail to tail. 

These observations help to understand the dynamic 
process of filament merging that leads to the persistent 
patterns in the open flows at intermediate Da. With the 
consideration of periodic boundary conditions, it is not 
difficult to understand also the process ending excitation 
in the closed flow at not too large Pe in terms of the 
annihilation of halves of two-humped filaments. 

What seems to escape from the picture is the process 
ending the persistent pattern in the open flows at large 
Da: because of the openness of the flow, there are always 
filaments that do not collide with others, but that receive 
the fresh reactants entering the system. In this case, as in 
the situations of very large Pe, it seems that strain fluctu- 
ations are essential to understand the process of deexcita- 
tion. We have run model (15) with A randomly changing 
in time and found that, for large enough fluctuations, 
nonvanishing correlation time, and values of Da in the 
two- hump regime: the induced width fluctuations even- 
tually end with the decay of the filament, in very much 
the same way as seen in the baker models simulations. 
We speculate that the larger width fluctuations origi- 
nated by filament collisions in situations such as those 
in Fig. 10 (right panel) would amplify still more the ef- 
fect of unsteady strain and help to eliminate excitation. 
Strain fluctuations in the baker model (periodic applica- 
tion of contraction followed by periods without strain) 
are an artifact of the discrete nature of the model. But 



in the two-dimensional simulations of Sect. Ill, and in 
real flows, strain fluctuations occur naturally and may be 
thus responsible for the excitation decay in the open flow 
at large Da. It is quite natural that this decay process 
only appears after the filaments develop the two-hump 
structure, since this implies the presence of refractory 
material. Nevertheless, a quantitative description of this 
process is still missing. 

V. CONCLUSIONS 

We have analyzed the behaviour of an excitable 
medium in the presence of open or closed chaotic flows. 
In both cases, three different regimes have been eluci- 
dated. The one at smaller Da, that is the dilution of the 
excitation at fast stirring, is analogous to what is found 
in the case of bistable chemical dynamics [13]. 

The most interesting regimes are found at intermedi- 
ate Da: in the closed flow coherent excitation of 
the whole system arises from the localized perturbation, 
whereas in the open flow the excitation remains indef- 
initely in the system. This last phenomenon was also 
found under bistable and in autocatalytic dynamics [13], 
as well as the excitation phase under the closed flow, that 
is, the growth of an excited filament that becomes space 
filling. What is distinct of the excitable dynamics is that 
excitation under the closed flow is a transient, so that the 
system finally recovers the rest state, at variance with the 
bistable and autocatalytic behaviour. It is striking that 
this recovery does not occur under the open flow in this 
intermediate Da range. 

Also a consequence of the recovery behaviour that 
characterizes the excitable dynamics, and that distin- 
guishes it from the otherwise rather similar bistable dy- 
namics, is the loss of coherence occurring at large Da. It 
manifests gradually under the closed flow, but as a sud- 
den disappearance of permanent excitation in the open 
case. 

A great part of our work has been devoted to the de- 
velopment of simplified models that help to understand 
the above regimes and the transitions among them. De- 
spite the strong approximations performed, these sim- 
ple models reproduce, at least qualitatively, the full two- 
dimensional numerical results. The first simplified model 
is based in the use of a baker map for the advection 
dynamics. It highlights the processes of stretching and 
folding as the basic flow mechanisms leading to the afore- 
mentioned chemical regimes. The baker model dynamics 
also suggest that transitions are linked to the proper- 
ties of individual filaments. Thus, an even simpler model 
is considered, where the stationary transverse profile of 
a filament is the main quantity under study. Most of 
the numerical observations can be understood within this 
framework, although some phenomena, specially those 
for which filament interactions seem relevant, would need 
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of more detailed modeling. 

Some geophysical observations have been already in- 
terpreted within the present framework [14]. It would 
be of great interest to perform experiments of chemical 
dynamics under well-controlled stirring to observe the 
different scenarios predicted here. 
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